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In this paper we discuss two issues, the inter-comparison of four mixed layer mesoscale parameteriza- 
tions and the search for the eddy induced velocity for an arbitrary tracer. It must be stressed that our 
analysis is limited to mixed layer mesoscales since we do not treat sub-mesoscales and small turbulent 
mixing. 

As for the first item, since three of the four parameterizations are expressed in terms of a stream func- 
tion and a residual flux of the RMT formalism (residual mean theory), while the fourth is expressed in 
terms of vertical and horizontal fluxes, we needed a formalism to connect the two formulations. The stan- 
dard RMT representation developed for the deep ocean cannot be extended to the mixed layer since its 
stream function does not vanish at the ocean’s surface. 

We develop a new RMT representation that satisfies the surface boundary condition. As for the general 
form of the eddy induced velocity for an arbitrary tracer, thus far, it has been assumed that there is only 
the one that originates from the curl of the stream function. This is because it was assumed that the tracer 
residual flux is purely diffusive. 

On the other hand, we show that in the case of an arbitrary tracer, the residual flux has also a skew 
component that gives rise to an additional bolus velocity. Therefore, instead of only one bolus velocity, 
there are now two, one coming from the curl of the stream function and other from the skew part of 
the residual flux. In the buoyancy case, only one bolus velocity contributes to the mean buoyancy equa- 
tion since the residual flux is indeed only diffusive. 

Published by Elsevier Ltd. 


1. Introduction 

In this work we discuss two issues: the intercomparison of four 
available mixed layer ML mesoscale parameterizations and 
whether the eddy induced velocity for buoyancy can also represent 
tracers other than buoyancy, for example, passive tracers such C0 2 , 
CFC, etc., that form part of climate studies. We study mixed layer 
mesoscales only, with no reference to sub-mesoscales and small 
scale turbulent mixing which require parameterizations not dis- 
cussed here. 

As for the first item, three of the four parameterizations are ex- 
pressed in terms of the stream function 4* and residual flux F, of 
the residual mean theory RMT formalism (Aiki et al„ 2004; Ferrari 
et al„ 2008, 2010, cited as A4, F8, 10), while the fourth one (Canuto 
et al„ submitted for publication, cited as Cl 1 ) is expressed in terms 
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of the vertical-horizontal mesoscale fluxes (F v , Fh ). 1 To carry out 
the model intercomparison, we need to translate the Cl 1 (F„, F H ) for- 
mulation into the corresponding formulation in terms of (*F, F r ). 


1 Killworth (2005, K5) was the first to argue that in the ML flows occur mostly 
within on horizontal planes and thus the natural representation of the mesoscale 
tracer flux is in terms of the horizontal F H and vertical F v component. I<5 solved the 
linear mesoscale dynamic equations and showed that F v is a skew flux, i.e., its 
divergence yields an horizontal advection with a bolus-like velocity u*. while F H is of 
the diffusion type with a mesoscale diffusivity ic M . Given the linear character of the 
model, 1<5 was unable to derive the strength of either u* and/or of k m . Recently, K5’s 
analysis was extended to include the non-linear terms (Cl 1 ) and the form of F H and F v 
for an arbitrary tracer in terms of the large scale fields was derived. When the tracer 
was the buoyancy field, the parameterization was assessed in several ways, e.g., z- 
profile of the eddy kinetic energy vs. WOCE data, surface eddy kinetic energy vs. T/P 
altimetry data, dependence of the vertical flux on the mean velocity against eddy 
resolving simulation data, etc. The ('P, F r ) representation has the advantage of 
facilitating the matching with the ocean interior at the bottom of the mixed layer 
while the (F v , F H ) representation has a different advantage. Since the dynamic 
equation for the EKE (eddy kinetic energy, see e.g., Boning and Budich, 1992) shows 
that F v acts as a source of EKE, one can model the surface eddy kinetic energy by 
averaging the vertical buoyancy flux F„ over the mixed layer and then assess the 
result against the T/P data (Scharffenberg and Stammer, 2010). 
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Since the standard form of (*F, F r ) cannot be used in the ML since it 
does not satisfy the boundary condition 4*(0) = 0, we developed a 
new RMT formulism valid in the ML. The final result, Eq. (21), ex- 
presses (F v . Fh) in terms of (*F, F r ). 

The results of the four models intercomparison can be summa- 
rized as follows: (a) the F8 bolus velocity does not entail ML re- 
stratification which is known to exist, (b) the A4, F10 bolus veloc- 
ities induce re-stratification but at the lowest order in the small- 
ness parameter fi/H, A4.F10 are not different (/i, H are the ML and 
ocean depths, respectively), (c) using the *P of A4, F8,10, we con- 
struct the corresponding F r ’s; we reproduce the F8 result while 
the F r of A4, F10 are new since they were not given in the original 
work, (d) in both A4, F10, F r is diffusive with a mesoscale diffusivity 
that vanishes at the bottom of the ML, as expected, (e) however, 
their values at the surface is /i/H times smaller than the value ob- 
tained by Zhurbas and Oh (2003), and finally, (f) only the Cl 1 mod- 
el accounts for wind and mean flow, which affect both the 
mesoscale fluxes and their kinetic energy. 

Concerning the parameterization of an arbitrary tracer, we ob- 
tain the following results. In addition to a diffusive component, 
the residual flux exhibits a new feature, a skew component, which 
gives rise to an additional bolus velocity. There are therefore two 
mesoscale advection terms: one due to the bolus velocity originat- 
ing from the stream function and the other from the bolus velocity 
originating from the skew part of the residual flux. The common 
assumption that there is only one bolus velocity is therefore no 
longer tenable (in the case of buoyancy, only the bolus velocity 
from the stream function contributes to the mean buoyancy 
equation). 

2. Inapplicability of the standard RMT to the ML 

Consider the model independent dynamical equation for the 
mean buoyancy b = -gpd'p (e.g., Ferreira et al„ 2005, Eq. (1)) 


Fr(fc) 


F(b) ■ Vb 
|Vb| 2 


Vb 


V H b + N 2 e z 
\Vb\ 2 


[Fh-Vh/j + FvN 2 ] 


( 4 ) 


where e z is the vertical unit vector and N is the Brunt-Vaisala fre- 
quency. The first term in (2) representing the isopycnal component 
of the buoyancy flux, has the form of a skew flux (Griffies, 1998) and 
its divergence leads to an advection: 

V • (<F x Vb) = U + ■ Vb (5a) 

where U + is the eddy induced or bolus velocity: 

U + = Vx>P, V H ■ u + + &W+ = 0 (5b) 


It follows that the mesoscale buoyancy flux F(b) contributes to Eq. 
(1) as an advection and a diffusion: 


V ■ F(b) = U + ■ Vb + V ■ F r 


(5c) 


In a fully adiabatic ocean, that is, one with no diabatic ML, the resid- 
ual flux F r is negligible and thus the main mesoscale effect is repre- 
sented by the first, advective, term in (5c). As McDougall and 
McIntosh (2001) showed, at the ocean’s surface the stream function 
satisfies the boundary condition 4*(0) = 0. We concentrate on the 
horizontal component of this condition: 

'Fh(O) = 0 (6a) 

since (6a) ensures the vanishing of the vertical component of the 
eddy induced velocity (5b) at the ocean’s surface: 

w+ (0) = e 2 ■ [V H x 'Ph(O)] = 0 (6b) 

Furthermore, condition (6a) leads to the vanishing of the vertical 
component of the isopycnal flux at the surface: 

'F(O) x Vb ■ e z = 0 (6c) 


Since the vertical component of the full flux also vanishes at the 
surface: 


<9 t b + U ■ Vb + V ■ F(b) = -V ■ F SM (b) - 9 Z F SS + Q. (1) 

Here, U = (u,w) is the mean velocity and F(b) = U'b' is the 3D 
mesosc ale buoyanc y fl ux with horizontal-vertical components 
F H (b) = u'b',F v (b) = w'b', U' = (u', w') is the mesoscale velocity field 
and V is the 3D nabla operator; the overbar stands for an ensemble 
average. 2 The first and second terms on the rhs represent the contri- 
bution due to sub-mesoscales and small scale (ss) turbulence (which 
we write for completeness but which we do not treat in this work), 
while Q stands for sources and sinks. The RMT decomposition of the 
buoyancy flux F(b) into isopycnal-diapycnal components is as fol- 
lows (Andrews and McIntyre, 1976; Treguier et al., 1997; Plumb 
and Ferrari, 2005; Ferreira et al., 2005; F8): 

F(b) = WF = ¥ x Vb + F r (b) (2) 

where the stream function pseudo-vector and the residual vector 
flux F r are defined as follows: 

▼ = - F( ^ |2 V5 = ^2 [( f vV H b - N 2 F„) x e z - F h x V H b] (3) 


2 The average in the mesoscale flux (1) stems from the non-linear term in the 
equation for the instantaneous buoyancy field when averaged over a coarse 

resolution grid cell of horizontal scale ~100km. However, averaging over the grid 
cell in (1 ) is not sufficient in, say, testing a mesoscale flux parameterization using high 

resolution simulations. The reason is that the diagnosed fluxes are random functions 
of the large scale fields which are the ones averaged over the grid cell. To obtain 
deterministic functions from high resolution data, one needs to ensemble average the 
random functions. If the large scale fields are stationary and/or homogeneous for 
sufficiently long time and/or within large area, the diagnosed instantaneous fluxes 
averaged over the grid cell may be further averaged over the corresponding time and / 
or space intervals instead of ensemble averaging. Below, we imply instantaneous 
averages over the grid cell together with ensemble averages. 


F(0) • e z = F v (0) = 0 (7) 

an analogous boundary condition must be satisfied by the residual 
flux: 


F r (0) • e z = 0 (8) 

How does the presence of the diabatic ML affect conditions (6)-(8)? To 
answer the question, we consider 4* H near surface. Since the last 
term in (3) does not contribute to 'F H and the first term is very 
small, strictly, it vanishes at z = 0 because of (7), we consider the 
term: 


N 2 (z)F H (z) x e z 


JV 4 (z) + |V H b| 2 


(9a) 


Near the surface, the horizontal flux F H does not vanish as it follows 
from the observational result by Zhurbas and Oh (2003): 

F h ( 0) = -K s V H b, K s = k m ( 0) = a/C ,/2 (0) (9b) 


who arrived at it using data from the Global Drifter Program/Surface 
Velocity Program. In (9b), k m (z) is the mesoscale diffusivity and 
C= 1.02 ±0.13; furthermore ^ = min(r d i R ) where r d is the Rossby 
deformation radius and L R is the Rhines scale. The surface mesoscale 
eddy kinetic energy I<( 0) can be obtained from the T/P data (Scharf- 
fenberg and Stammer, 2010). Since in the ocean, even near the sur- 
face, one has s < 1 (typical isopycnal slopes below the ML are of the 
order of 1CT 3 while in the ML they are about an order of magnitude 
larger) from (9a, b) we conclude that condition (6a) and therefore 
(6b, c), are not satisfied. In addition, the interpretation of the diver- 
gence of the skew flux (5a) as an advection becomes problematic. In 
fact, in the limit scl, from (9a, b) we obtain that the vertical bolus 
velocity becomes: 
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w+(0) = e z ■ [V H x ’Fh(O)] = -V H • (K s s) (9c) 

which is a diffusive not advective term. 

Next, let us consider the limiting case N 2 ( 0) = 0. Although from 
(9a) it follows that the (6a) is satisfied, the following difficulty arises. 
In the ML, the variables |V H b| ~ 1CT 8 s~ 2 and F H (z) = -K M V H b. are 
almost z-independent while N 2 (z) increases from its assumed sur- 
face value JV 2 (0) = 0 to N 2 (-h) ~ 10~ 5 s~ 2 * at the bottom of the ML. 
It follows that the function ¥ H vs. N 2 given in relation (9a) has a max- 
imum at N 2 = | V H b|, at a depth denoted by z m: 

z = z m N 2 = V H b ~ 1CT 8 s~ 2 < N 2 (-h) ~ 1CT 5 s~ 2 (10a) 

Substituting this relation into (9b, a), we obtain: 

'Fh(Ziti) ~ - x e z (10b) 

1 |V H b| 

Using this result, we compute the average of the bolus velocity from 
the surface to the depth z m : 

(u + ) z = e z x < dzV > Zm = e z x ¥(z m )/z m = — k s z 1 (10c) 

1 |V H b| 

Thus, near the surface, the bolus velocity is estimated to be: 

z = 0: (u + ) Zm = 1 ^ ~ 10 ms -1 (10d) 

z z m 

where we have taken k s ~ 10 3 m 2 s~’and z m ~ 100 m. Eq. (lOd) is 
the difficulty we alluded to earlier since it represents an unphysical 
result. In the next section we show that the following modification 
of the RMT functions: 

¥,F r — > ¥,F r (lOe) 

solves the problem by ensuring the compliance of conditions (6a, c) 
and how the stream function ¥ now varies smoothly from zero at 
the surface to |¥(-h)| ~ 1 m 2 s _1 at the bottom of the ML. This 
behavior dispenses with the need to use artificial tapering schemes 
(Griffies et al., 2005) whose arbitrariness resulted in a widely vary- 
ing characteristics of the simulated flows (Gnanadesikan et al., 
2007; F8). It must also be noted that over the years, several authors 
argued for the need of the new (tilde) variables, e.g., Treguier et al. 
(1997), Held and Schneider (1999), Plumb and Ferrari (2005) and F8 
even though the tilde itself was not explicitly used. 

3. Formulation of an RMT for the mixed layer 

In order to satisfy ¥(0) = 0 while keeping (2), we carry out the 
following transformations: 

¥^¥ = ¥ + A. F r ^F r = F r -AxVS (11) 

Here, A is an arbitrary function of space and time for whose deter- 
mination we follow the following principles. Since there is no rea- 
son to modify the RMT formulated for the deep ocean, we search 
for a ¥ that is almost identical to ¥ in the deep ocean but satisfies 
the condition ¥(0) = 0. This in turn requires that the function A 
vanish in the deep ocean and that it cancels ¥ at the surface. Since 
the solution of the problem is not unique, additional conditions 
must be imposed. We choose ¥ to be a 2D pseudo-vector and fol- 
lowing F10, we introduce the vector Y: 

¥=Yxe z (12a) 

Thus, the eddy induced velocity defined in Eq. (5b) with ¥ -> ¥, is 
expressed as follows: 

u + = <9 Z Y, w+ = -VY (12b) 

The general form of A that satisfies both the above requirements is 
the following: 


|V H b| 2 A = N 2 F r xe z + F H x V H b + D (13) 

where D is an arbitrary 2D pseudo-vector that satisfies the condi- 
tion D(0) = 0 and that vanishes in the interior. It is worth mention- 
ing that F8 used only the first term in (13). Next, we consider the 
D = 0 case and show that it does not satisfy physical requirements 
that, in addition to the boundary conditions, must be met. We then 
suggest an expression for D that solves the problem. 

3.1. The D = 0 case 

Substituting (11) and (13) into Eqs. (3) and (4), we obtain the 
following results 3 : 

Y = -|V6r 2 N 2 F‘ H r + |V H br 2 FvV H h (14) 

F r = Ff, + | V H br 2 FvN 2 V H E + |V H 5| 2 |Vbr 2 F[; (15) 

At the surface, both F v and F[) vanish, the former by construction, 
the latter because relation (9b) shows that at z = 0 the flux F H is di- 
rected along the horizontal gradient of the mean buoyancy and 
thus, has no transverse component. From relations (12a) and (14), 
it then follows the desired relation ¥(0) = 0 which, in terms of 
Y(z), becomes Y(0, -~H) = 0. Next, if we multiply (14) scalarly by 
V H b, we obtain the following two relations: 

F v = Y'-V H b. N 2 F f H r = -|Vb| 2 Y tr (16) 

which, once substituted into (15), yield the following result: 

F r = Fh + N 2 Y f + |V H b| 2 hT 2 Y tr (17) 

This relation connects the stream function to the residual flux and 
shows that the parameterization of the two variables cannot be 
chosen independently. 4 In the next sections we will discuss the 
use of analogous relations to study the A4, F8, 10, Cl 1 models. Final- 
ly, we note that in relations (16) and (17), when N 2 vanishes, there is 
a singularity unless we have that: 

Y tr (N 2 = 0) = 0 (18) 

As we discuss in Section 4, this condition is indeed satisfied in F8 
but not in A4 and F10. Thus, to avoid singularities in the latter cases, 
we need to include an appropriate function D into (13). 

3.2. The D ^ 0 case 

Let us consider the choice: 

|Vb| 2 D = J\r 2 |V H b| 2 e z x f£ (19) 

which is allowed only if F|J vanishes not only at z = 0 for any N but 
also when N = 0 for any z. In particular, in the Cll model discussed 
in Section 4, this condition is indeed satisfied. Then, instead of rela- 


tions (14) and (15), we now have: 

Y = -N- 2 F t H r +F v |V H hr 2 VHb 

(20a) 

F r = F' H +F v JV 2 |V H 5r 2 V H 5 

(20b) 


Thus, instead of (16) and (17), we now have the final relations: 

¥ = (Y tr + Y f ) x e z , Y tr = -FT 2 F[[, 

Y l = ^%S7 H b, F r = F' H (b) + N 2 Y‘ (21) 


3 To simplify the notation, we have introduced the transverse and longitudinal 

components of an arbitrary horizontal vector V which are defined as fol- 
lows:V = V tr + V‘, |v H b| V' = (V V H b)V H b. 

4 It is instructive to consider a zonal flow when Vnb = d v be v , u'b' = v'b'ey. Then, 

from (14) and (15) we obtain H* = F r = F r e y , F r = v'b' + N 2, F. The last relation 

links the stream function to the residual flux and coincides with the results of Held 

and Schneider (1999) who suggested it using heuristic arguments. 
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The problem of finding the stream function and the residual flux corre- 
sponding to a given vertical and horizontal fluxes, is solved. Specifi- 
cally, the transverse component of F H and F v yield the stream 
function ¥ ; the longitudinal component of F H then yields the resid- 
ual flux. Relations (21 ) show that the stream function and the resid- 
ual flux depend on each other and so do their parameterizations. To 
highlight the difference between *F and ¥, we rewrite them as 
follows: 

Old RMT : ¥ = aN- 4 F v V H b x e z - ai\T 2 F H x e z - ahT 4 F H x V H b 

( 22 ) 


ML: u + *-h- C Mp4 (26) 

2 hNj 

which does not depend on z. Using N 2 = db/dz, from Eq. (1) we 
obtain: 

^=-d z (u+.Vfa)-d z (V-F r ) + --- (27) 

Because both u + and Vb are z- independent, the eddy induced advec- 
tion does not contribute to the ML stratification. This result was one of 
the motivations for the F10 model we discuss next. 


New RMT : ¥ = aN- 4 F v V H b x e z - J\T 2 F H x e z - al\r 4 (F H • s) 

x Vh b x e z (23) 

where a = (1 + s 2 ) -1 and a = s~ 2 . Within a down-gradient model the 
last term in (22) vanishes. Since F v ( 0) = 0, the first term in (22) van- 
ishes while the second term does not vanish since F H (0) # 0, see 
relation (9b). Thus, ¥(0) ^ 0, as already discussed. On the other 
hand, using a down-gradient model of the type F H (b) = -fc M V H h. 
the last two terms in (23) cancel each other out while the first term 
vanishes at z = 0. This results in ¥(0) = 0, as required. 


4. Comparison of four parameterizations 

4.1. F8. model 


This model assumes that between the ML of thickness h and the 
deep interior, there is a transition layer (TL) of thickness D < h. In 
the ML, the stream function ¥ is assumed to be linear in z and to 
vanish at the surface. In the TL, ¥ is parameterized as a quadratic 
function of z. At the boundary between ML and TL, ¥ and its z 
derivative are required to be continuous while at -Z\ = H + D , 
which is the boundary between TL and the interior, ¥ and its z- 
derivative are required to match those of the GM model; the hor- 
izontal buoyancy flux is assumed to be of the down-gradient type 
in all three regimes ML, TL and deep ocean. This assumption al- 
lowed the authors to get rid of the vertical and transverse compo- 
nents of ¥ before the transformation (11) was carried out. This 
feature, in turn made the second and third term in (13) redundant 
and thus the definitions of the modified stream function and resid- 
ual flux have the forms given by Eqs. (12a), (14), and (15). On this 
basis, the model yields (see Eqs. (25) and (28) of F8): 


Y(z) = (CgmG(Z) ^ , F r = — K GM 

Nf 


1 - G(z) 


N 2 (z) 

Nf 


V H b 


(24) 


where G(z) is given in Eq. (26) of F8. As one can see, (24) satisfies 
relation (17) in which Y tr = 0, as we discussed below Eq. (17). To 
make (24) in the ML and TL more transparent, we write the function 
G(z) in the limit D<gh which is in accordance with F8’s estimates of 
D and h given in their Fig. 5a and b. In addition, we assume that 
within the TL, the z-derivative of N 2 is constant (in the ML, F8 as- 
sume that N 2 <c N 2 is constant). This implies that in the TL, 
D _1 = -N/ 2 d z N 2 . Adopting a rigid lid approximation, relation (26) 
of F8 becomes: 


-h^z^O: G(z)f«-^^; -(h + D)<z<-h: G(z) 

3 z (z + h) 2 
2 h 2D 2 


(25) 


Substituting (24) and (25) into (12b) and since the main contribu- 
tion to the z-derivative comes from the differentiation of G(z), we 
obtain that in the ML: 


4.2. A4 and F10 models 

These models suggest the following differential equations for 
Y(z): 

A4 : C-’4Y = -V H fa, F10 : (c 2 d 2 zz -N 2 )Y = -K CM V H b 

(28a, b) 

Since at any depth, the solutions of (28a, b) depend on V H b at all 
depths and since V H E has different directions at different depths, 

Y has both longitudinal and transverse components. As it follows 
from (17), in the D = 0 case the transverse component of Y leads 
to singularities in F r . Thus, the A4, F10 models need the choice 
(19) which avoids the singularities, as Eqs. (21) show. 

As one can see, the difference between the two models is the term 
N 2 (z). In the opinion of the authors of F10, in the A4 model without 
such a term, ‘‘information about the background stratification, and 
hence about the mixed layer, would be lost." Below, we study analyti- 
cally the A4, FI 0 models and reach the following conclusions: (1) to 
the first order in h/H, in both A4 and FI 0, the ML stream function does 
not depend on the stratification N 2 and, (2) such a dependence is 
present in F10 but it comes into play only to higher orders in h/H 
(h, FI are the ML and ocean depths). Thus, FI 0 differs from A4 in those 
locations where h ~ FI, while in the rest of the ocean the two models 
are quite similar. To carry out the analysis, we employ 
JV 2 (z) = N 2 = const, in the ocean’s interior and JVj(, L < N 2 . As for the 
horizontal buoyancy gradients in the ML and interior, we assume 
that they are different Vh^m # V H b i but z-independent in their 
respective regimes. F10 suggest that the speed c in (28b) be chosen 
as the first baroclinic velocity and thus, to the main order in h/H, 
we have c = N t Flln. We have found that the ML solution of Eq. (28b) 
with the boundary condition Y(0, - H ) = 0 is given by 5 : 

Y = -KgmC(oA + BC), C = nz/H (29) 

where the two vectors A, B and the constant a are defined as 
follows: 

N 2 A = V H hi. JV 2 B = ^Vh£>m, a = tanh n/2 (30) 

In (29), we kept the term since it is the only one that contributes 
to the re-stratification of the ML, as we show below. In the A4 mod- 
el, a = 7t/2. In both models, a was determined by solving the bound- 
ary-value problem between the surface and the bottom. To the main 
order in h/H, these results do not depend on the ML stratification N. The 
bolus velocity and its z-derivative contributing to the ML re-strati- 
fication, are obtained from (12a) and (29): 

- N 2 

u + = -7iK'G M H _1 (aA + 2BQ, -9 z u + ■ V H b = 4n 2 \B\ 2 K CM -j > 0 

H 

(31) 


5 In the Eady case corresponding to N 2 (z) = Nf = const. throughout the whole 
ocean treated in F10, the solution of Eq. (28b) is given by Eqs. (38) and (39) of F10. 
When expanded in power of z/H, the result Y Eady « -K/ GM AC(tanh 7r/2 + C/2) coincides 
with (29) since in the Eady case 2A = B. 
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In contrast to the F8 model, the second of (31) leads to a ML re- 
stratification. 

Though in A4 and FI 0 the residual flux was not presented, it can 
be found on the basis of our last relation (21) and the solutions of 
(28a, b). Though A4 and F10 do not give the residual flux, it can be 
found using their stream function. First, we study F[, and consider 
the top of the ocean interior where F r « F r « 0. From the fourth of 
(21), we then obtain: 


(c) the residual flux was not presented in A4 and F10. We have 
shown that in both A4.F10 it is of the diffusive type with a 
mesoscale diffusivity that vanishes at the bottom of the ML, 

(d) however, the value of the surface diffusivity in such models 
is roughly h/H times smaller than what Zhurbas and Oh 
(2003) have derived using drifter data. 

4.3. Cl 1 model 6 


z = zi F' H «-N 2 Y f 


(32a) The Cl 1 results for the horizontal and vertical buoyancy fluxes 

are: 


On the other hand, since inside the ML the variables u', b' and V H b 
are almost z-independent, in the region z Ss i\ we may use (9b) and 
extend it all the way to the bottom of the ML. Thus, using the 
approximate relation: 

z 3= z, : F„ « -k s V h 5 (32b) 


and substituting it into (32a), we can estimate the surface value of 
the mesoscale diffusivity: 


K s « N 2 Y f • V H b/ V H b 


(32c) 


where all the variables are computed at the bottom of the mixed 
layer. Substituting (32b, c) into the last of (21) we obtain: 


A4.F10: F r = -fc r (z)V H f> 


(33a) 


K r 


( Y-Vh b |V H b,| 2 N 2 (zA 

\ Y |V H fa| 2 N? J 


K s 


z N 2 (z)\ 

h^r) 


(33b) 


where in the last step we used Eqs. (29) and (30). Thus, the residual 
flux Eq. (33) is represented by a horizontal diffusion with a z-depen- 
dent mesoscale diffusivity which equals k s at the surface and van- 
ishes at z = -h. To evaluate k s , we substitute (29) and (30) into 
(32c) with the results: 


Ks(F10) 

TCgm 


7t[tanh(7t/2)](h/H), ~ U 2 (h/H) 

Kgm ^ 


(34a) 


Since relations (34a) were obtained for the particular case when in 
the interior the buoyancy gradient is z independent, they should be 
viewed as estimates which we indicate by These results show 
that in the case of a sufficiently shallow ML where /j/H<tc 0.2, we 
have: 

A4, F10 : K s < K C m (34b) 

While the commonly used value Kgm ~ 10 3 m 2 s _1 (e.g., Ferreira 
et al., 2005) is consistent with the second relation (9b), use of the 
typical values r d = 30 km and K = 1CT 2 m 2 s~ 2 in relation (34b) en- 
tails a surface diffusivity much lower than what is observed. This 
conclusion deserves some comments. We arrived at (34a, b) thanks 
to relations (21) that allows us to relate the stream function of the 
A4-F10 models to the horizontal flux which in turn we related to 
(9b) derived from observations. Since at the time A4 and F10 were 
presented, relation (21) was not known, the translation of 4* into a 
horizontal flux and then the comparison with the data based rela- 
tion (9b), was not possible, conclusions (34) could not be reached 
and went unnoticed. In conclusion: 

(a) the F8 eddy induced velocity does not entail ML re-stratifica- 
tion which is known to exist, 

(b) the A4 and F10 eddy induced velocities yield re-stratification 
but, in spite of their formal difference, at the lowest order in 
the small parameter h/H, neither model exhibits a depen- 
dence of the re-stratification on the ML N 2 (z). Therefore, 
A4 and F10 are not different, 


F H (b) = -K M V H b, F v (b) = -k V H b (35a) 

where: 

k = k m zF(z) (35b) 

/r d F(z) = (u + u)xe z , u(z) = u(z)-u d 

zu(z) = f u (z')dz', (•) ee f /C 1/2 (z)dz / f /( 1/2 (z)dz (35c) 

Jo J-H / J-H 

The variable u d represents the mesoscale “drift velocity” which in 
the f-plane has the following expression: 

u d = (u) - i/r d e z x < d z s > (35d) 

The above relations exhibit an interesting physical feature. Since 
mesoscale eddies move with their own drift velocity u d , the effec- 
tive mean velocity is the one in the frame co-moving with the mes- 
oscales and that is why the mean velocity enters as in the second 
relation in (35c). Finally, the mesoscale diffusivity is given by: 


Km (z) =/(u,/C)fK I/2 (z) 


(36a) 

where: 



/(u,fC)= l+^|u(z)| 2 

, f =min( r d ,L R ) 

(36b) 


When the eddy kinetic energy K is larger than that of the mean flow, 
see Fig. 7 of Scharffenberg and Stammer (2010), the function f(u,K) 
is close to unity. 

To present the Cl 1 parameterization in terms of the RMT for- 
malism, we first remark that from (35) and (21) we have: 

F„ = 0. Y tr = 0, Y f = -k 1 = - K ' VH f V H b (37) 

|V H b| 2 

where k‘ is the longitudinal component of k, see footnote 3. Using 
(35b), Y acquires the form: 

Y = -zk M F ' V , V H b (38a) 

|V H b| 2 

To link this result with the CM model, we rewrite it as follows: 

Y = T(z)Y GM (-h), T(z) = ■ V H b (38b) 

|VhO| 

where the function T(z) satisfies the conditions 

7(0) = 0, T(-h) = 1 (38c) 


6 In the case of buoyancy, the Cll parameterization was assessed in several ways: 
z-profile of the eddy kinetic energy vs. WOCE data, surface eddy kinetic energy vs. T/P 
altimetry data, dependence of the vertical flux on the full mean velocity field (not only 
its geostrophic part) which is affected by wind, overall assessment against eddy 
resolving simulation data, etc. 
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While the first relation is obvious, the second condition comes from 
the fact that at the bottom of the ML, the smallness of the measured 
diapycnal diffusivities (Ledwell et al„ 1998, 2011) implies the van- 
ishing of the diapycnal flux: 

F A =Fv + N- 2 F H V H b (39) 

Two considerations are in order here: (1) relations (38b, c) are anal- 
ogous to the default tapering scheme (Griffies et al., 2005; Gnanad- 
esikan et al., 2007) with the difference that the tapering function T(z) 
is no longer arbitrary but given by the model, as the second relation 
in (38b) shows and, (2) we did not impose that at the bottom of the 
ML the stream function Y matches that of GM, and yet it does, as the 
first of (38b) shows. As for the residual flux, use of (35a, b) in the last 
of (21), yields the following results: 

F r = -K r V H b, K r = jc m [1 -T(z)], T(z) = -^—j-T(z) (40a) 

N (-h) 

in which the diffusivity K r vanishes at the bottom of the ML. Eqs. 
(38b) also show that in the Cll model, the profile of k t and the 
tapering function T(z) relate to one another while within the default 
tapering scheme they are independent. Finally, we notice that in 
most of the ML, N 2 is sufficiently small and thus from (40a) we 
have: 

z > -h K r k' m (40b) 

a relation that we use in the next section. 

For further analysis, it is convenient to decompose F into F 5 and 
F 2 components that are governed by the deep ocean and the ML 
respectively: 

F = F, + F 2 , F, = -rVud x e z , F 2 (z) =/“ V(u + u) x e 2 

(41a) 

and consequently: 

Y = Y,+Y 2 , Yi >2 = —KmZF\ 2 (41b) 

Since both terms in (35d) are of the same order of magnitude 
A, where A is given in (30), use of the relation r d = N]H/7r[/l 
yields (the subscript I refers to the interior): 

Yi~-K M J(n.A)n, n=|V H b|^V H 5 (41c) 

Thus, the first term of (29) and (41b) are rather similar. As for Y 2 , we 
distinguish three cases. 

4.3.3. No wind, geostrophic flow 

Using the thermal wind relation, 0 = zf ’e z x V H bivi = 2u and 
r d = we obtain the following expression (the suffix “g” 

stands for geostrophic): 


where u gs is the geostrophic component of the mean surface 
velocity. 

4.3.3. Ekman layer 

The presence of an Ekman layer, represented here for simplicity 
by an Ekman spiral, contributes a term given by: 

Y e = |V H br 2 F E V H b, F e = A(j/)t • V H b + fi(*7)V H 5 x t • e z (42b) 

where t is the wind stress, i/ = z/c5 E , <5 E is the thickness of the Ekman 
layer and: 

A(tj) = (1 - e'' cos r\) + r/e' 1 (sin i] - cos rf) 

’ * p ° (42c) 

B(tj) = (T /|/|) [» 7 (cos r\ + sin rf) +sin;/]e'' 

Thus, the complete form of the stream function is given by: 

Y = Y 1 +Y 2g + Ysv+YE (43) 

New terms 

While the terms Yi,Y 2g in (41c,d) are similar to the corresponding 
terms in A4-F10, the terms Y sv , Y E given in (42a, b) are absent in 
A4, F10. One should distinguish Y sv , Y E from the terms appearing 
in the mean buoyancy equation due to the presence of wind depen- 
dent terms in the mean velocity. While the latter stem from the di- 
rect effect of the wind stresses on the mean momentum equation, 
Y sv , Ye originate from the non-linear interaction of mesoscales and 
mean velocities in the dynamic equations for the mesoscale 
velocities. 

In summary, the main difference between Cll and A4, F10 is 
that in the former case, the surface diffusivity coincides with zc M , 
Eq. (36), in agreement with observations, whereas in A4, F10 the 
surface diffusivity is much smaller than what the data indicate. 
The Cll model yields profiles of T, F r that are analogous to those 
of the default tapering schemes with the difference that the pro- 
files are no longer arbitrary but determined within the model. 

5. The case of an arbitrary tracer 

In the previous sections we considered the ML mesoscale 
parameterization for buoyancy. However, OGCMs time step tracers 
such as temperature, salinity, C0 2 , CFC, etc. which cannot be repre- 
sented as a buoyancy field. In what follows we show that the eddy 
induced velocity for buoyancy and arbitrary tracers are different. 
We begin with the equation for an arbitrary tracer: 

dtZ + U- Vt + V-F(t) = — V • Fsm — <9 Z F SS + Q (44a) 

where the rhs has the same meaning as in Eq.(l ). Using the relations 
(35) for the case of a tracer: 

Fh(t) = -K M V H T, Fv(t) = -k • V H T (44b) 


Y - 3 k VHfaM r 2 
*2g - ~ 9 . ,2 e > 

^ iVj 


: = 


7TZ 

JT 


we have: 

(41 d) 

F(t) = -k m V h t - (k ■ V H t)e z 


(44c) 


which, after using (30), is close to the second term in (29). Thus, in Next, consider the second term in the rhs of (44c) which we rewrite 

the absence of wind stresses, A4, F10 and CIO yield similar results. as follows: 


4.3.2. Wind 

In the presence of wind, the situation changes significantly 
since there are wind stresses and horizontal gradients of the sur- 
face pressure which, in turn, result in a geostrophic surface velocity 
which, as shown in Cl 1 , dominates the production of ML eddy ki- 
netic energy. In the presence of a geostrophic surface velocity (de- 
noted by the subscript sv), the term Y 2 has an additional term that 
is not present in the A4, F10 models: 


-(k • V H T)e z = -(k • Vi)e z = — (k x e z ) x Vi - k — 
_ , . dr 

— Fskew(^) ^ Qz 

where we have defined the skew flux: 

Fskew(T) = -(K X e z ) X Vt 
Eq. (44c) then becomes: 


„ Km . 

Ysv = — 2 2 x ‘ 

2fr d 


(42a) F(r) = -K m V h T - k— + F skew (T) = F diff (T) + F skew (T) 


(45a) 

(45b) 


(45c) 
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Substituting (45c) into (44a), we obtain: 

d t X + (D + U„) ■ Vt + V • Fdiff (T) = - V ■ Fsm - <9 Z F SS + Q. (45d) 
where the eddy induced velocity U, =(u»,w*) is given by: 

dK _ . . _ . 

u * = _ az' w * = Vh K ( 45e ) 

We must remark that in the buoyancy case, Eq.(45d) has the form: 

dtb + (D + U + ) • Vb + V ■ Fdiff(b) = -V • F SM - <9 Z F SS + Q. (46a) 
where (V • U + = 0): 

, <9Y dK 1 . f . 

u + = — , w + = V H • k 1 (46b) 

Since: 


u 


dK 1 _ 

dz ~ 

it follows that 

u, = u + + u« 


dK 

dz 


dK tt 

dz 


u, - u. 


(46c) 


(46d) 


Therefore, the so-called “residual velocities" for buoyancy and arbi- 
trary tracers in Eqs. (46a) and (45d) are given by: 


U r es ( b) = U + U + , U res ( T) = U res (b) + U 


(46e) 


The common assumption (e.g., Ferreira and Marshall, 2006) that 
there is only one bolus velocity u + is no longer valid since in the case 
of an arbitrary tracer the true bolus velocity is (46d) in which the 
new component u„ contributes only to tracers other than buoyancy 
since in the latter case: 


u« ■ V H b = 0 (46f) 

It is important to note that the same results are obtained using 
the RTM formalism. In fact, substituting Eqs. (45c, b), (12a) and (37) 
into the definition of the tracer residual flux: 

F r (t) = F(r) - >P x Vt (47a) 

we obtain: 

F r (T) = Fduf(T) + F skew (T) (47b) 

where: 

Fskew(T) = -K tr X Vt (47 c) 


Therefore, within the Cl 1 model, the residual flux has not only a dif- 
fusive component as all previous models, but also a skew one which 
gives rise to the additional bolus velocity U„. Substituting 
F(t) = T x Vt+ F r (T) into (44a), together with (12a), (37) and 
(47b, c), we obtain the equation: 

<9rT + (U + U + U„) • Vt + Vh ■ Fdiff = — V • Fsm — <9 Z F SS + G 

(47d) 

which, because of (46d), is identical to (45d). 

6. Summary and conclusions 

In this work, we have compared four recent ML mesoscale 
parameterizations. In Cll, the parameterization was presented 
in terms of the horizontal and vertical mesoscale fluxes F H and 
F v . While F h is given by a down-gradient diffusion, F v is a skew 
flux whose divergence yields an advection. The other three 
parameterizations, F8,10 and A4 were formulated in terms of a 
stream function and a residual flux of the residual mean theory 
RMT, though the last two modeled only the stream function. For 


this reason, in the first part, we compared the four parameteriza- 
tions for the buoyancy flux only. To compare the Cll results 
with the other three models, we had to translate the F h and F v 
fluxes into the mean residual theory formalism, i.e., a stream 
function and a residual flux. We began with discussing the fact 
that the definition of the latter used for the adiabatic ocean inte- 
rior, must be modified in the diabatic ML in order to satisfy the 
homogeneous boundary condition of the stream function. The 
appropriate modifications and the new stream function and 
residual flux were discussed in Section 3. The expressions for 
the modified ’F and F r in terms of the V-H fluxes are given in 
Eq. (21). These relations allowed us not only to compare Cll 
with F8,10 and A4, but will be instrumental in matching the 
ML mesoscale parameterization with the one in the adiabatic 
ocean. We shall treat this problem in subsequent work. 

In the absence of the wind stresses, we have shown that the 
stream functions of A4, F10 and Cll are similar. In the case of 
strong winds, the Cll model yields a stream function that in- 
cludes the Ekman and surface velocities whereas these two con- 
tributions are absent in the A4, F8, 10 models. The inclusion of 
the Ekman flow is important since eddy resolving simulations 
(Maltrud et al„ 1998) showed that the surface eddy kinetic en- 
ergy is contributed quite substantially by the Ekman flow. As 
for the re-stratification effects of mesoscales, we have shown 
that to the main order in the parameter h/H (which is small 
everywhere except in deep convective regimes), the models A4, 
F10 give the same result and neither of them exhibits a depen- 
dence on the ML stratification. 

An interesting result emerges from the new stream function 
and residual flux, namely that they are related to each other, see 
Eq. (21), and, thus, parameterization of one affects the other, a fea- 
ture that did not exist in the deep ocean where the residual flux is 
negligible. 

We have checked that relation (21 ) is satisfied in the F8 and Cl 1 
models while the A4, F10 do not provide an explicit parameteriza- 
tion of the residual flux. We have shown that in order to satisfy 
(21 ), the A4 and F10 models require a surface mesoscale diffusivity 
that is much smaller than the results of the observations by Zhur- 
bas and Oh (2003). 

Concerning the parameterization of an arbitrary tracer, we have 
found the following: 

(a) for OGCMs that employ vertical-horizontal fluxes in the ML, 
the appropriate arbitrary tracer equation is given by an 
equation of the form (44a) with the fluxes given by Eqs. 
(44b). 

(b) for OGCMs that employ the RMT in the ML, the arbitrary tra- 
cer equation is Eq. (45d) 

(c) within the TRM formalism, and in the case of an arbitrary 
tracer, the eddy induced velocity u + computed as the curl 
of the stream function does not represent the full eddy 
induced velocity which is given by u„. 
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